#include <cmath>
#include <cstdio>

#include "geo2d.h"

namespace geo2d {

    double Distance(Point2d a, Point2d b) {
        return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
    }
    double Distance2(Point2d a, Point2d b) {
        return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
    }
    void getTriangleCircumcenter(Point2d a, Point2d b, Point2d c, Point2d& center, double& radius) {
        double x1 = a.x, y1 = a.y;
        double x2 = b.x, y2 = b.y;
        double x3 = c.x, y3 = c.y;
        double A = x1 - x2;
        double B = y1 - y2;
        double C = x1 - x3;
        double D = y1 - y3;
        double E = (x1 * x1 - x2 * x2 + y1 * y1 - y2 * y2) / 2;
        double F = (x1 * x1 - x3 * x3 + y1 * y1 - y3 * y3) / 2;
        double det = B * C - A * D;
        center.x = -(D * E - B * F) / det;
        center.y = -(A * F - C * E) / det;
        radius = sqrt((center.x - x1) * (center.x - x1) + (center.y - y1) * (center.y - y1));

        double radius1 = sqrt((center.x - x2) * (center.x - x2) + (center.y - y2) * (center.y - y2));
        double radius2 = sqrt((center.x - x3) * (center.x - x3) + (center.y - y3) * (center.y - y3));
        if(fabs(radius - radius1) > 1e-6 || fabs(radius - radius2) > 1e-6) {
            printf("radius mismatch radius: %lf, radius1: %lf, radius2: %lf\n", radius, radius1, radius2);
        }
    }
    
    double crossProduct(Point2d a, Point2d b) {
        return a.x * b.y - a.y * b.x;
    }
    
    // 计算向量 AB 与向量 AC 的叉积
    // 返回值 >0 表示 C 在 AB 的左侧，<0 表示 C 在 AB 的右侧，=0 表示三点共线
    double crossProduct(const Point2d &a, const Point2d &b, const Point2d &c) {
        return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
    }

    // 判断点 q 是否在线段 pr 上（包含端点），适当考虑误差 EPS
    bool onSegment(const Point2d &p, const Point2d &q, const Point2d &r) {
        if (q.x <= std::max(p.x, r.x) + EPS && q.x >= std::min(p.x, r.x) - EPS &&
            q.y <= std::max(p.y, r.y) + EPS && q.y >= std::min(p.y, r.y) - EPS)
            return true;
        return false;
    }

    // 计算三个点的相对方向
    // 返回 0：共线，1：逆时针，2：顺时针
    int orientation(const Point2d &p, const Point2d &q, const Point2d &r) {
        double val = crossProduct(q - p, r - p);
        if (fabs(val) < EPS)
            return 0; // 共线
        return (val > 0) ? 1 : 2;
    }

    // 判断两条线段是否相交
    bool segmentsIntersect(const Point2d &p1, const Point2d &p2,
                           const Point2d &p3, const Point2d &p4) {
        // 分别计算四个方向
        int o1 = orientation(p1, p2, p3);
        int o2 = orientation(p1, p2, p4);
        int o3 = orientation(p3, p4, p1);
        int o4 = orientation(p3, p4, p2);
        
        // 一般情况：若 p1p2 与 p3p4 的两个端点的方向不同，则相交
        if (o1 != o2 && o3 != o4)
            return true;
        
        // 特殊情况：若共线则判断是否有端点在另一条线段上
        if (o1 == 0 && onSegment(p1, p3, p2)) return true;
        if (o2 == 0 && onSegment(p1, p4, p2)) return true;
        if (o3 == 0 && onSegment(p3, p1, p4)) return true;
        if (o4 == 0 && onSegment(p3, p2, p4)) return true;
        
        return false;
    }

    // 利用参数化方程求直线交点（适用于非平行情况）
    // 直线1：A -> B，直线2：C -> D
    // 若两直线有唯一交点，则返回 true，交点保存在 intersect 中；否则返回 false
    bool lineIntersection(const Point2d &A, const Point2d &B,
                          const Point2d &C, const Point2d &D,
                          Point2d &intersect) {
        Point2d r = B - A;
        Point2d s = D - C;
        double rxs = crossProduct(r, s);
        if (fabs(rxs) < EPS) {
            // 直线平行或共线，无法求唯一交点
            return false;
        }
        double t = crossProduct(C - A, s) / rxs;
        intersect = A + r * t;
        return true;
    }

    bool InCircle(Point2d p1, Point2d p2, Point2d other) {
        Point2d center = Point2d((p1.x + p2.x) / 2, (p1.y + p2.y) / 2);
        double radius2 = Distance2(p1, center);
        double dist2 = Distance2(other, center);
        return dist2 <= radius2;
    }
    
    double triangleArea(Point2d a, Point2d b, Point2d c) {
        return 0.5 * fabs(crossProduct(Point2d(b.x - a.x, b.y - a.y), Point2d(c.x - a.x, c.y - a.y)));
    }

    bool isPointInsideTriangle(Point2d a, Point2d b, Point2d c, Point2d p) {
        // Compute area of full triangle ABC
        double A = triangleArea(a, b, c);

        // Compute areas of sub-triangles PAB, PBC, and PCA
        double A1 = triangleArea(p, b, c);
        double A2 = triangleArea(a, p, c);
        double A3 = triangleArea(a, b, p);

        // If the sum of A1, A2, and A3 equals A, then P is inside the triangle
        return fabs(A - (A1 + A2 + A3)) < 1e-6;        
    }

    bool isPointInTriangle(const Point2d &a, const Point2d &b,
                           const Point2d &c, const Point2d &pt) {
        // 利用叉积法：如果 pt 与三角形各边组成的向量叉积均同号，则在三角形内
        double cp1 = crossProduct(a, b, pt);
        double cp2 = crossProduct(b, c, pt);
        double cp3 = crossProduct(c, a, pt);
        bool hasNeg = (cp1 < -EPS) || (cp2 < -EPS) || (cp3 < -EPS);
        bool hasPos = (cp1 >  EPS) || (cp2 >  EPS) || (cp3 >  EPS);
        return !(hasNeg && hasPos);
    }

    std::pair<bool, int> checkQuadrilateralConvexity(const Point2d &p1, const Point2d &p2,
                                                    const Point2d &p3, const Point2d &p4) {
        // 计算4个连续三点的叉积
        double cross1 = crossProduct(p1, p2, p4);
        double cross2 = crossProduct(p2, p4, p3);
        double cross3 = crossProduct(p4, p3, p1);
        double cross4 = crossProduct(p3, p1, p2);
        // printf("cross1: %f, cross2: %f, cross3: %f, cross4: %f\n", cross1, cross2, cross3, cross4);
        
        constexpr double eps = 1e-6; // cant set 1e-9 because of zero value check
        bool allNonNegative = true, allNonPositive = true;
        double arr[4] = { cross1, cross2, cross3, cross4 };
        for (int i = 0; i < 4; i++) {
            if (arr[i] < -eps) allNonNegative = false;
            if (arr[i] > eps)  allNonPositive = false;
        }
        
        // printf("allNonNegative: %d, allNonPositive: %d\n", allNonNegative, allNonPositive);

        if (allNonNegative || allNonPositive) {
            // 四个叉积同号（允许零值），说明四边形凸
            return {true, -1};
        } else {
            // 若存在正负不同，则为凹四边形
            // 统计正负个数
            int posCount = 0, negCount = 0;
            for (int i = 0; i < 4; i++) {
                if (arr[i] > eps)
                    posCount++;
                else if (arr[i] < -eps)
                    negCount++;
            }
            // printf("cross1: %f, cross2: %f, cross3: %f, cross4: %f\n", cross1, cross2, cross3, cross4);
            // 以多数符号为准
            double majoritySign = (posCount >= negCount) ? 1.0 : -1.0;
            // 对于每个叉积，不符合多数符号的，认为其对应的中间点为凹点
            // 顺序：cross1 -> p1, cross2 -> p2, cross3 -> p4, cross4 -> p3
            if (cross1 * majoritySign < -eps) {
                return {false, 1}; // p1为凹角
            }
            if (cross2 * majoritySign < -eps) {
                return {false, 2}; // p2为凹角
            }
            if (cross3 * majoritySign < -eps) {
                return {false, 4}; // p4为凹角
            }
            if (cross4 * majoritySign < -eps) {
                return {false, 3}; // p3为凹角
            }
            // 理论上不会到此处
            return {true, -1};
        }
    }

    // 平面内四边形查找
    void searchPointOnSurface(const std::vector<Point2d> &points, const std::vector<int> &is_convex, const Point2d &p, int &x, int &y, int nx, int ny) {
        // 逐一检查每个四边形
        if (points.size() != (nx + 1) * (ny + 1) || is_convex.size() != nx * ny) {
            printf("size mismatch\n");
            return;
        }

        for (int j = 0; j < ny; j++) {
            for (int i = 0; i < nx; i++) {
                // 取当前单元4个角点（按顺序）
                // example: 2 - - - 3
                //          |       |
                //          |       |
                //          0 - - - 1
                Point2d p1 = points[j * (nx + 1) + i];
                Point2d p2 = points[j * (nx + 1) + i + 1];
                Point2d p3 = points[(j + 1) * (nx + 1) + i];
                Point2d p4 = points[(j + 1) * (nx + 1) + i + 1];
                // printf("p1: (%f, %f), p2: (%f, %f), p3: (%f, %f), p4: (%f, %f)\n", p1.x, p1.y, p2.x, p2.y, p3.x, p3.y, p4.x, p4.y);

                // is_convex: 0 - convex, 1-4 - concave point index
                if (is_convex[j * nx + i] == 0) {
                    // 凸四边形
                    // printf("convex branch\n");
                    double cross1 = crossProduct(p1, p2, p);
                    double cross2 = crossProduct(p2, p4, p);
                    double cross3 = crossProduct(p4, p3, p);
                    double cross4 = crossProduct(p3, p1, p);
                    // printf("cross1: %f, cross2: %f, cross3: %f, cross4: %f\n", cross1, cross2, cross3, cross4);
                    if ((cross1 < 0 && cross2 < 0 && cross3 < 0 && cross4 < 0) ||
                        (cross1 > 0 && cross2 > 0 && cross3 > 0 && cross4 > 0)) {
                        x = i;
                        y = j;
                        return;
                    }
                } else {
                    // 凹四边形
                    // printf("concave branch\n");
                    Point2d arr[4] = { p1, p2, p3, p4 };
                    int concaveIndex = is_convex[j * nx + i] - 1; // 取出凹角顶点索引（0~3）
                    int oppIndex   = 3 - concaveIndex; // 对应凹角对面的顶点
                    int nextIndex  = concaveIndex ^ 1; // 对应凹角的下一个顶点
                    int prevIndex  = concaveIndex ^ 2; // 对应凹角的前一个顶点

                    // 分别构成两个三角形：
                    // 三角形1： (arr[concaveIndex], arr[nextIndex], arr[oppIndex])
                    // 三角形2： (arr[concaveIndex], arr[oppIndex], arr[prevIndex])

                    double tri11 = crossProduct(arr[concaveIndex], arr[nextIndex], p);
                    double tri12 = crossProduct(arr[nextIndex], arr[oppIndex], p);
                    double tri13 = crossProduct(arr[oppIndex], arr[concaveIndex], p); // new line

                    double tri21 = crossProduct(arr[concaveIndex], arr[oppIndex], p); // new line
                    double tri22 = crossProduct(arr[oppIndex], arr[prevIndex], p);
                    double tri23 = crossProduct(arr[prevIndex], arr[concaveIndex], p);

                    bool inTri1 = (tri11 < 0 && tri12 < 0 && tri13 < 0) || (tri11 > 0 && tri12 > 0 && tri13 > 0);
                    bool inTri2 = (tri21 < 0 && tri22 < 0 && tri23 < 0) || (tri21 > 0 && tri22 > 0 && tri23 > 0);

                    if (inTri1 || inTri2) {
                        x = i;
                        y = j;
                        return;
                    }
                }
            }
        }
        printf("point not found\n");
    }
}; // namespace geo2d
